Cramér-Rao Bound Optimized Subspace Reconstruction in Quantitative MRI

We extend the traditional framework for estimating subspace bases that maximize the preserved signal energy to additionally preserve the Cramér-Rao bound (CRB) of the biophysical parameters and, ultimately, improve accuracy and precision in the quantitative maps. To this end, we introduce an approximate compressed CRB based on orthogonalized versions of the signal's derivatives with respect to the model parameters. This approximation permits singular value decomposition (SVD)-based minimization of both the CRB and signal losses during compression. Compared to the traditional SVD approach, the proposed method better preserves the CRB across all biophysical parameters with negligible cost to the preserved signal energy, leading to reduced bias and variance of the parameter estimates in simulation. In vivo, improved accuracy and precision are observed in two quantitative neuroimaging applications, permitting the use of smaller basis sizes in subspace reconstruction and offering significant computational savings.

A common approach for reconstructing such undersampled dynamic data is to use different sub-sampling patterns for each contrast and utilize redundant information between contrasts, e.g., with low-rank subspace reconstruction techniques [5], [10]- [16].This approach assumes that the different contrasts lie in some low-dimensional subspace and reconstructs coefficient images corresponding to the basis functions.This subspace is usually precomputed by taking the singular value decomposition (SVD) of a simulated dictionary of signals over the expected range of tissue parameters, e.g., using the Bloch equations or the extended phase graph method [2], [17]- [19].
In qMRI, high-quality image reconstruction is desirable, but not in itself as-outside of nonlinear inversion approaches [8], [20]- [23]-it is only an intermediate step to obtaining quantitative parameter maps [24], [25].The most common approach is to fit the biophysical model directly to the coefficient images, e.g., with dictionary matching [4], [12], [26], non-linear least squares, kernel regression [27], or neural-network-based [28]- [30] fitting.In this paradigm, the subspace coefficients of each voxel are the measurements used for parameter estimation.Consequently, the variance of an unbiased estimator is bounded by the Cramér-Rao lower bound [31] (CRB) of the measured coefficients rather than the uncompressed data (e.g., a time series of images).In this article, we analyze the conservation of the CRB during the projection of the signal onto the subspace and present a strategy for optimizing the subspace to mitigate the associated increase in CRB.
CRB analysis has a rich history in MRI where it has often been used to design the optimal experimental parameters, such as the flip angle train or the undersampling pattern [32]- [40].In the context of subspace reconstruction, recent work analyzed the g-factor for a CRB-optimized pulse sequence and fixed basis to identify the optimal subspace size for reconstruction [41].In this work, we optimize the basis itself, hypothesizing that, while minimum variance unbiased estimation is generally unachievable in qMRI, improved CRB preservation leads to improved conditioning of the parameter estimation process.While reduced CRBs might be expected simply to improve parameter precision, we show that, in some situations, it can also lead to improved parameter accuracy.This is consistent with viewing the CRB as a measure of the local curvature of the log-likelihood function over parameter space [24], where higher curvature can reduce bias [42].In joint estimation problems, lower CRBs for difficult-to-estimate parameters can also decrease the bias sensitivity for easier-toestimate parameters due to reduced information coupling [43].
In this work, we propose a method for optimizing linear bases to preserve the CRB of the biophysical parameters in addition to the signal energy by leveraging the CRB's geometric interpretation [44].We incorporate orthogonalized versions of the signal's derivatives into a basis optimization scheme as an approximation for the compressed domain CRB, which we analyze in detail.We demonstrate our approach's ability to improve the accuracy and precision of the parameter maps with smaller subspace sizes in silico and in vivo in two qMRI neuroimaging applications: a) a two-pool qMT model [45]- [47] using a hybrid-state sequence [48], and b) singlecompartment T 1 and T 2 mapping using the MRF inversion recovery fast imaging with steady-state precession (MRF-FISP) technique [26].This is an expansion of our previous work [49], where we considered only the qMT model and did not analyze the approximation made to the compressed CRB.

A. Problem Description
In this article, we adopt the following linear subspace model, which does not consider the MR imaging process: where s ∈ C N T is the measured signal with N T time points or data frames, U ∈ C N T ×Nc is the subspace of rank N c where N c ≪ N T , ′ is the conjugate transpose, c ∈ C Nc is the coefficient vector representing the contribution of each basis element to the spin evolution, and ϵ ∼ N (0, σ 2 I) is white Gaussian noise (I is the identity matrix).Note that , the low-rank operator does not change the noise covariance if U encompasses an orthonormal basis.
As the coefficients c are the measurements received by the parameter estimator, we use this model to generate an expression for the compressed CRB in Section II-D.

B. Traditional SVD Basis
The subspace for a qMRI experiment can be estimated a priori.Using an appropriate signal model, N s simulated signal evolutions (or fingerprints) for the expected range of parameters can be stacked column-wise to form a dictionary matrix S ∈ C N T ×Ns .Then, U is the solution to the following sample principal components analysis problem: where UU ′ represents the projection onto ⟨U⟩ (angular brackets ⟨•⟩ denoting the column span) and ∥ • ∥ F denotes the Frobenius norm which captures the signal energy loss between the compressed and original fingerprints.Eq. ( 1) is efficiently solved by taking the first N c left singular vectors of the SVD of S to form U [50], which we refer to as the "traditional SVD" approach.
C. The Cram ér-Rao Bound and its Geometric Interpretation For an unbiased estimator, the CRB is the minimum variance of a parameter's estimate.For Gaussian noise, each parameter's CRB is proportional to the corresponding diagonal entry of the inverse Fisher Information Matrix F −1 = (J ′ J) −1 , where J is the Jacobian matrix whose columns are the signal's derivatives with respect to the parameters θ, i.e., j i ≜ ∂s/∂θ i where we consider, for the moment, the uncompressed signal s as the measurements.
For a one-parameter model, the CRB is simply proportional to the inverse ℓ 2 -norm of the derivative: As described in Ref. [44], this simple notation can be translated to the multi-parametric case by considering the angle between the different signal derivatives.Specifically, the CRB in the multi-parametric case is inversely proportional to the ℓ 2norm of the derivative w.r.t.θ i after removing all components parallel to the linear space spanned by the derivatives w.r.t.all other model parameters, which we denote using the matrix J i .Concretely, the uncompressed Cramér-Rao bound B u (θ i ) of the parameter θ i can be written as where P Ji ≜ J i (J ′ i J i ) −1 J ′ i denotes the projection matrix onto ⟨J i ⟩.To highlight the similarity to the CRB of the singleparameter model, we defined the orthogonalized derivative j i,⊥ as the projection of the derivative j i onto the space orthogonal to ⟨J i ⟩, as shown in Fig. 1.For clarity we denoted the projection matrix twice, noting that they are symmetric and idempotent, i.e., (I − P Ji ) ′ (I − P Ji ) = I − P Ji .Eq. ( 2) enables the calculation of a parameter's CRB from the ℓ 2norm of its orthogonalized derivative alone-without requiring a matrix inverse-which helps formulate a minimization procedure.It also reveals that the orientations and norms of the spaces ⟨j i ⟩ and ⟨J i ⟩ are what determine the CRB.Increasing the dimension of ⟨J i ⟩, e.g. by adding additional model parameters, usually decreases ∥j i,⊥ ∥ 2 and consequently increases the CRB.Decreasing the dimensionality, e.g. by fixing certain parameters, usually has the opposite effect.

D. The Cram ér-Rao Bound for Subspace Coefficients as Measurements
In this section, we modify Eq. (2) to incorporate subspace modeling.In this case, we consider the subspace coefficients c as the measurements used for parameter estimation.The modification consists of simply multiplying the signal derivatives with U ′ , yielding the exact compressed CRB (B ec ) where the latter equality follows from definition of the projection above.Eq. ( 4) shows that the CRB is increased relative to Eq. ( 2) if the projection onto ⟨U⟩ reduces the angle between ⟨j i ⟩ and ⟨J i ⟩.However, the CRB remains unchanged if ⟨J⟩ ∈ ⟨U⟩, i.e., if the basis captures the span of the entire Jacobian matrix.While this condition is fulfilled if the signal S ∈ ⟨U⟩, this is generally unachievable when N c ≪ N T , meaning there are also no guarantees surrounding the preservation of the Jacobian.In fact, maximizing the preserved signal energy in Eq. (1) equates to optimizing only for the signal's derivative with respect to the scaling M 0 without accounting for any other derivative and their geometric relations.

E. Cram ér-Rao Bound Optimized Bases
To both represent the signal evolution accurately and preserve the CRB of the biophysical parameters, we propose to optimize a linear basis Û by incorporating both objectives in a composite cost function, given by minimize ) where J ⊥ ∈ C N T ×NsNp is the matrix whose columns are the N p orthogonalized derivatives of interest for each dictionary fingerprint, each normalized to have unit energy.This formulation permits optimizing for only a subset of the model parameters (i.e., the parameters of interest for any given application) while still considering a fit of the full model in computing the CRB.Normalization of each orthogonalized derivative allows us to measure the proportion of the CRB lost in compression and weights each parameter in the cost equally, where it would otherwise be dominated by the most difficult-to-estimate parameters.While we have not done so here, we note that it is still possible to incorporate weightings for the different parameters into Eq.( 5).λ ∈ [0, 1] controls the convex combination of the signal energy and CRB losses in the overall cost.Defining λ in this way enables us to probe the behavior of the bases over a fixed range of λ values across different pulse sequences, which may have different signal or derivative amplitudes.Setting λ = 0 reduces Eq. ( 5) to the traditional SVD approach in Eq. ( 1).
The second term in Eq. ( 5) captures the approximate CRB loss (∆B ac ) and can be rewritten as where the right-hand term is the ratio of the uncompressed CRB to an approximate compressed CRB (B ac ), defined as Comparing Eq. ( 4) and Eq.(7) shows that the approximation lies in the projection onto ⟨J i ⟩ instead of ⟨UU ′ J i ⟩.The exact compressed CRB loss (∆B ec ) analogous to Eq. ( 6) is In general, the following relation holds: where the latter inequality follows directly from because the projection onto ⟨U⟩, a subspace of C N T , can only decrease the angle between ⟨J i ⟩ and j i .If the angle is zero after compression, i.e. ⟨UU ′ J i ⟩ = ⟨UU ′ j i ⟩, B ac would be finite while B ec is infinite.The following relationship also follows directly from Eq. ( 9): meaning that minimizing the approximate CRB loss ∆B ac does not in general guarantee minimization of the exact CRB loss ∆B ec .However, this is less likely to occur for sufficiently large values of λ as ⟨J⟩ is increasingly well-preserved within the CRB-SVD basis, and hence B ac ≈ B ec .
While we would ideally like to directly minimize ∆B ec in Eq. ( 8), we introduce ∆B ac to utilize the established and numerically preferable SVD framework to optimize the basis functions.We first compute the orthogonalized derivatives {j i,⊥ } NsNp i=1 , which do not depend on U. Thereafter, Eq. ( 5) can be solved by simply combining the two loss terms and performing an SVD on the horizontally concatenated matrix where ŨNc denotes the first N c columns of Ũ. Eq. (12) shows that the primary difference between the CRB-SVD and the traditional SVD approach is the explicit inclusion of the model parameter's orthogonalized signal derivativesrather than only the signal (i.e., ∂s/∂M 0 )-in the subspace estimation process.This helps preserve the orientations of the subspaces ⟨j i ⟩ and ⟨J i ⟩ ∀i, and therefore the CRB, at some cost to the preserved signal energy.
The SVD of D naturally yields optimal (within the approximation), orthogonal bases which need only be computed once per value of λ and can be retrospectively truncated according to the choice of N c .Note that because the column dimension of D scales linearly with N s N p , assuming N T < N s N p , the memory complexity O (N T N s N p ) and compute time complexity O N 2 T N s N p of the CRB-SVD scale linearly with N s N p [51].

A. Pulse Sequences
The first qMRI application we consider is the hybridstate sequence described in Ref. [48] designed to extract the parameters of a 2-pool qMT model [45]- [47], i.e. a complexvalued scaling M 0 , the fractional semi-solid spin-pool size m s 0 , the free spin-pool relaxation rates R f 1 , R f 2 , the exchange rate R x , the semi-solid spin-pool relaxation rates/times R s 1 , T s 2 , and the field inhomogeneities B 0 and B + 1 (9 total parameters).This hybrid-state sequence is optimized to minimize the CRB for and T s 2 in individual 4s long cycles with antiperiodic boundary conditions [9], [48].We use 3D radial koosh-ball k-space sampling with a 2D golden means pattern [52] reshuffled to minimize eddy current artifacts [53].

B. Data Simulation
For the qMT model, we simulated a dictionary of approximately 600,000 fingerprints with the generalized Bloch framework [47] using three Gaussian distributions representing brain tissue (WM+GM), fat, and cerebrospinal fluid (CSF), as outlined in Table I. 67% and 33% of the total fingerprints were used for basis calculation and testing, respectively.
For MRF-FISP, we used Bloch simulations to compute a dictionary of 281,250 fingerprints over the Cartesian grid shown in Tab.I, accounting for spoiling across the slice profile [56] by taking the complex average of 1324 isochromats (5300 for cerebrospinal fluid) [57].

C. Basis Optimization
For the qMT model, we consider 6 parameters of interest (N p = 6): 1 , and T s 2 , and treat M 0 , B 0 , and B + 1 as nuisance parameters which are only used indirectly in constructing J ⊥ to calculate the orthogonalized derivatives for the parameters of interest.For MRF-FISP, we consider T 1 and T 2 (N p = 2) while treating M 0 as a nuisance parameter.The a We use three Gaussian distributions, denoted by N (m, s) with mean m and standard deviation s, corresponding to typical values in grey and white matter, fat, and cerebrospinal fluid at 3T [26], [48], [54], [55].The scripts denote truncation limits and U denotes a uniform distribution.The percentile brackets denote the relative size of each tissue type.b MRF-FISP values are simulated on a discrete grid (min:stepsize:max).
signal derivatives for a given fingerprint are orthogonalized with respect to one another using QR factorization.
For the hybrid-state pulse sequence analyzed here, a sign change of B 0 entails simply a complex conjugation of the signal.Consequently, we expect real-valued bases with a symmetric B 0 distribution.To avoid approximation errors when randomly drawing samples, we complement each fingerprint with its complex conjugate to ensure that this symmetry is not broken and that we obtain real-valued bases.
For both applications, we perform the CRB-SVD using 10 evenly spaced values λ ∈ [0, 1), omitting λ = 1 due to poor signal fidelity leading to severe image artifacts.SVDs were performed with 20GB memory distributed over 10 CPU threads (Intel Skylake 6148, Santa Clara, California, USA) for MRF-FISP and 700GB memory/30 CPU threads for the qMT model on our institution's computational cluster.The latter took an average of 7 hours for each value of λ.
To evaluate the CRB-SVD method, we compare B ac (Eq.( 7)) to B ec (Eq.( 4)) as a function of N c and λ.We also define the ratio as a proxy for the optimality of the CRB-SVD basis.While the CRB-SVD computes a globally optimal solution for preserving B ac , this solution is close to globally optimal for preserving B ec only when R ≈ 1.Let us denote the CRB-SVD solution with the values B 0 ac and B 0 ec , and the values for a global minimizer of B ec with B 1 ec and B 1 ac .Then, since ec provides a tight bound on B 1 ec , the exact compressed CRB we would ideally like to minimize.Note R ≤ 1 by Eq. ( 9).

D. Simulation Experiments and Parameter Fitting
To evaluate the impact of the CRB-SVD basis on parameter fitting, for both our test applications described in Section III-A we evaluate the bias and standard deviation of white matter parameter estimates across 1000 noisy measurements.For the qMT model we used the values m s 0 = 0.2, R f 1 = 0.52/s, R f 2 = 12.9/s, R x = 16.5/s,R s 1 = 2.97/s, and T s 2 = 12.4µs [48], while for MRF-FISP we used T 1 = 810ms and T 2 = 25ms.The signal was simulated using the same methods described in Section III-B and compressed using the traditional SVD and CRB-SVD bases.We added Gaussian noise with an assumed SNR ≜ |M 0 |/σ = 50.
For the qMT model, we used the non-linear least squares (NLLS) estimator using the Levenberg-Marquardt algorithm [58] initialized with the ground-truth values.For MRF-FISP, we used dictionary matching over the same Cartesian grid described in Section III-B, where the stepsize was chosen to be smaller than the expected noise in the parameters to better approximate maximum likelihood estimation.Dictionary elements were compressed to the same subspace and matching was performed directly on the compressed coefficients [12].

E. In Vivo Imaging Experiments
For the qMT sequence, we scanned the whole brain of a healthy subject on a 3T Biograph mMR (Siemens, Erlangen, Germany) with a 32-channel head coil and a 256mm isotropic FOV.To facilitate reconstruction with large N c values-which are used for validation and have comparably high memory demands-we used a lower-than-usual spatial resolution of 1.6mm isotropic resolution and 30 cycles of the hybrid-state sequence (similar to a multishot acquisition) for 12.6 min scan time.For the MRF-FISP sequence, we acquired a single slice through another healthy volunteer's brain on a 3T Prisma MRI scanner (Siemens, Erlangen, Germany).We used a 32-channel head coil, an FOV=256mm x 256mm, voxel size 1x1x4mm, and 3 cycles (3 radial spokes per frame), for approximately 1 min of scan time.Informed consent was obtained for both subjects in agreement with our IRB's requirements.
For both sequences, we reconstructed coefficient images for various combinations of {λ, N c } in Julia [59] following the low-rank inversion approach [12], [14].We used a Toeplitz approximation of the non-uniform FFT [60], [61] for computational efficiency with non-Cartesian data and sensitivity encoding using coil maps calculated with ESPIRiT [62].For the qMT sequence, we used a fixed locally low-rank penalty [63], [64] strength and 250 iterations of FISTA [65] to suppress artifacts and noise, followed by voxel-wise NLLS fitting with a maximum of 500 iterations.For MRF-FISP we used the conjugate gradient algorithm followed by voxel-wise dictionary matching as described in Section III-D.

A. Quantitative Magnetization Transfer
Fig. 2 shows that increasing CRB-weightings (λ) leads to an improvement in CRB preservation at some cost to the preserved signal energy.This improvement is more pronounced for small N c : from λ = 0 to λ = 0.5 for N c = 15, the average ∆B ec decreases from 0.81 to 0.34-equating to a 72% decrease in the average CRB across all parameters-which notably is better than the preserved CRB for {λ = 0, N c = 30}.λ ∈ [0.4,0.6] provides nearly the maximal improvement in the CRB with minimal cost to the preserved energy for all N c .We observe that the R values (Eq.( 13)) achieved by the CRB-SVD bases have a cap of 0.85 for N c = 30-the limit of a computationally achievable reconstruction.This indicates the practical infeasibility of capturing the entire span of the signal derivatives within a compact basis for this application but also suggests good practical utility for the CRB-SVD basis: the traditional SVD basis for N c = 15 has R = 0.39, demonstrating its suboptimality in preserving B ec .
While B ac never quite approaches B ec , Fig. 3 investigates this approximation by plotting a regression across all test fingerprints.Note that all points lie below the reference identity line, consistent with Eq. ( 9).B ac correlates poorly with B ec for N c = 10-which is close to the number of parameters in this model-where deviations from the identity line are generally more severe for large CRB values.The correlation improves with increasing λ, and even more so with increasing N cwhich also significantly reduces the CRB, such that all test  points are captured within the plotted axis limits.
Fig. 4 shows that the bias and variance of NLLS-based m s 0 and R x estimates are reduced with the CRB-SVD (λ = 0.5) basis compared to the traditional SVD (λ = 0) for all N c .The difference is particularly pronounced for N c < 18, where the B ec reference shows that the estimation problem becomes particularly difficult with the traditional SVD basis (blue line) and NLLS deviates significantly from a minimum variance unbiased estimator.In particular, the variance is reduced relative to B ec at the cost of introducing significant bias.As there is no model mismatch in this simulation, this bias reflects a failure of the estimator rather than the biophysical model.However, both bias and variance remain stable for the CRB-SVD basis across all N c and both parameters.
In vivo, we similarly observe improved performance consistency using the CRB-SVD basis across all N c as shown in Fig. 5, where N c = 30 is used as a gold standard.The magnifications demonstrate the improved precision in m s 0 , most prominently for small N c .In contrast to the improved precision in m s 0 , the maps of the generally worse-conditioned R x parameter reveal mostly improvements in accuracy.The increasing CRB of R x for N c ≤ 18 results in a significant bias with the traditional SVD basis which is not seen with the use of the CRB-SVD basis.Here, increasing N c is of limited benefit, and increasing λ has a much greater impact on accuracy.Since R x is a comparably difficult-to-estimate parameter, the accuracy and precision of its estimate tend to guide the selection of the appropriate N c for this application.

B. MRF-FISP
Fig. 6 compares in vivo fits of R 1 derived from the traditional SVD and CRB-SVD bases with dictionary matching.Here, we observe that the CRB-SVD yields similar results for N c ≥ 4, with a benefit seen mostly at the extreme of N c = 3 where the traditional SVD's estimates are substantially biased.

V. DISCUSSION
Our results suggest that performing the SVD on the signals alone is suboptimal in preserving the CRB of the model's parameters, which depends on preserving the orthogonal components of the signal's derivatives.In this work, we proposed a method for incorporating CRB preservation into an SVD-based calculation of the basis.This approach utilizes a geometric interpretation of the CRB [44] and relies on an approximation of the compressed CRB.In general, we find that a stronger CRB weighting (increased λ) results in improved CRB preservation with a small cost to signal fidelity.Additionally, we find that a stronger CRB weighting also promotes a better approximation of the exact compressed CRB (B ec ) by the approximate compressed CRB (B ac ), improving confidence in the quality of the solution.In simulation, we found that when the number of measurements (N c ) was similar to the number of model parameters, the proposed method reduced both bias and variance of the parameter estimates, which will be further discussed in Section V-C.In vivo, the CRB-SVD basis showed some improvements at minimal N c when analyzing the MRF-FISP sequence and non-inferiority otherwise.For the qMT sequence, the proposed approach offered improved accuracy and precision at constant N c , or the ability to reduce N c while retaining accuracy and precision.
While we focused here on adding a CRB preservation term to the traditional SVD basis optimization objective, we emphasize that signal fidelity is still of fundamental importance.In the image reconstruction process, unrepresented signal energy causes inconsistency with the measured data, leading to artifacts scattered throughout the coefficient images that can bias the parameter estimates.Still, we have shown in this article that there is a range of λ-values where the improved  CRB outweighs the small cost in signal fidelity.This suggests that the slight increase in artifacts is less important than the improvements to the parameter estimation optimization landscape gained from reduced CRBs.

A. Practical Considerations
The proposed CRB-SVD method involves a simple modification to the traditional SVD basis calculation, requiring only precomputation and orthogonalization of the signal's derivatives with respect to the model parameters.These can be calculated analytically for many qMRI sequences, but automatic differentiation [37], [66] or physics-inspired neural networks [67] can also be used.The appropriate λ is application dependent, where a slightly higher λ was found to be suitable for the qMT application because the improved CRB, especially at low subspace sizes, outweighed the loss in preserved signal.We postulate that the selection of λ may, in general, also be sensitive to the signal variety-e.g. the variety of anatomy and tissue types-that needs to be accurately represented in the basis.Regardless, a moderate λ value of 0.3-0.5 generally appears to be a good starting point for other applications.Once computed, the CRB-SVD bases can be directly utilized with standard reconstruction pipelines, such as in BART [68] for typical subspace reconstruction tasks or even multitasking applications [5], [69], [70] where the basis for some dimension (e.g., T 1 relaxation) can be predetermined.
The practical benefit for subspace reconstructions in terms of computational and memory demands is significant, which can be limiting particularly for non-Cartesian 4D+ applications.Many iterative reconstruction algorithms require the computation and storage of the normal operators, whose size scales quadratically with N c .For example, in the qMT case, reducing 30 coefficients to 15 or even 10 reduces the memory requirements by a factor of 4 and 9, respectively.In practice, for a 160x150x120 matrix size (which corresponds to the full-brain coverage with 1.6mm isotropic resolution we used in this article)-assuming complex-valued singleprecision floating points and 2x oversampling in each dimension for a Toeplitz approximation-is the difference between 39GB/17GB and 154GB memory, which is often unavailable on a typical workstation.For 1.0mm isotropic resolution [48], the reduction is from 675GB to 169GB or 75GB memory, conferring feasibility on most high-performance computational platforms.A smaller number of coefficients also directly speeds up the calculation of the Toeplitz kernels and shortens the per-iteration computation time (the normal operators involve O(N 2 c ) complex-valued floating point multiplications per k-space data point), reducing the overall reconstruction time (assuming a fixed number of iterations).

B. Limitations and Drawbacks
While general superiority or even non-inferiority of the CRB-SVD basis is difficult to prove and is dependent on the figure of merit, we analyzed two extreme cases here to provide information across a wide spectrum of potential applications.We analyzed the MRF-FISP sequence as it is known for its very compact representation [12].For this reason, we did not expect major improvements and Fig. 6 largely confirms this expectation.At the other end of the spectrum, we analyzed a qMT sequence that has a less compact representation and showed that the CRB-SVD has substantial benefits in such a case.Many practical applications-and consequently the expected benefit-likely lie between these two illustrative examples.Nonetheless, we have identified other extreme cases from the literature: McGivney et.al. [12] used N c = 200 to represent the signal of the original MRF approach [4], and for such cases, we expect a substantial benefit when using the CRB-SVD over its traditional counterpart.
While the performance and optimality of the proposed CRB-SVD are straightforward to verify (c.f.Fig. 2), an important limitation is that we did not perform direct minimization of the exact compressed CRB.This would improve the optimality of the CRB-SVD bases with potential further improvements in representational efficiency (i.e., improved accuracy and precision at smaller subspace sizes).Direct minimization of the exact compressed CRB is not straightforward as the orthogonal projection onto ⟨J i ⟩ depends on U (Eq. ( 4)), and we need to additionally fulfill the orthonormal basis constraint (U ′ U = I) because our CRB expression assumes that the basis coefficients have no cross-talk, and bases that are not normalized can artificially decrease the CRB.This is a non-convex minimization over the complex Grassmannian likely requiring an iterative procedure, e.g.stochastic gradient descent, which will be the subject of future work.
In this article, we assumed a white Gaussian noise model and did not consider deviations under experimental conditions, e.g.due to parallel imaging [71] or compressed sensing.We note that the CRB is somewhat limited as a design criterion because it considers only local properties of the log-likelihood function and an optimistic CRB can sometimes mask poor model identifiability [72]; i.e., that the data can be explained by multiple combinations of parameters [42].For example, a model with poor identifiability would have many local minima which might cause convergence issues in gradientbased estimators [73]-which would have to be mitigated by the selection of appropriate constraints in parameter space.A Barankin-type bound, while more difficult to compute, would better capture the global structure of the statistical model and more accurately reflect realizable system performance in low SNR regimes [74] relevant in many qMRI applications, and is another avenue for future work.
One drawback of the proposed method is the memory requirements for performing the CRB-SVD, which as described in Section II-E scale linearly with the number of included orthogonalized derivatives and fingerprints.The sample size for the 9-parameter qMT model was dictated by the practical limitation of the maximum memory available on our computational cluster, though there was virtually no difference in the traditional SVD-derived basis with 5x the sample size.Additionally, for most applications, the basis functions have to be calculated only once, while the memory requirements for the reconstruction have to be fulfilled for every dataset.In future work, we will investigate more memory-efficient algorithms for performing large-scale SVDs [75].

C. Cram ér-Rao Bound: Precision or Accuracy?
The CRB is most often associated with the variance or noise of a parameter estimate.However, maximum likelihood estimators (e.g., the commonly used dictionary matching) are only asymptotically unbiased (or √ N c -consistent), i.e., for many measurements [76], [77].Least squares estimators (e.g., NLLS) are asymptotically unbiased only for white Gaussian noise models [76], [78].In vivo, where acquisition time is a major constraint, only limited measurements can be made.In qMRI, the number of measurements has often been similar to the number of model parameters (N c ≈ N p in the context of subspace reconstruction) tracing back to the original DESPOT method [1], which results in non-negligible bias.Unmodeled biophysical effects (e.g.unmodeled tissue compartments or diffusion in our qMT model [79]), regularization in the image reconstruction [80], and image artifacts all inexorably introduce further errors into the parameter estimation process.As a consequence, we practically cannot expect unbiased parameter estimation and the performance predicted by the CRB is rarely obtained.Nonetheless, the CRB is commonly used as a proxy for the "SNR-efficiency" or "conditioning" of the qMRI system [34], [39], and our work adopts this heuristic.Our simulations and in vivo experiments confirm that dictionary matching and NLLS estimates of the parameters are biased when N c ≈ N p , and particularly so with the traditional SVD basis, limiting the linear compressibility of advanced quantitative models that can practically be achieved.However, our results show that the use of the CRB-SVD basis improves both precision and accuracy in subspace reconstruction.
The observations about precision and accuracy with respect to the theoretical CRB made in this article are not new in the qMRI literature, and many papers report significant differences in both the mean and standard deviation of estimated parameters for CRB-optimized sequences [33]- [38].Notably, Zhao et.al. observe in simulated fully-sampled data an unbiased estimate of the relaxation times with the expected reduction in noise for a CRB-optimized sequence (note subspace reconstruction was not used in their work) [34].However, they find in the (spirally) undersampled regime that the parameter estimates are biased and the CRB-optimized sequence helps to improve both the accuracy and precision of the parameter maps.Unfortunately, limited conclusions about the relationship between the CRB and observed bias can be drawn from other studies involving CRB-based optimization of the scan parameters that show in vivo results [33], [35]- [38] due to bias introduced by unmodeled biophysical effects such as magnetization transfer.

VI. CONCLUSION
We present a method to incorporate CRB preservation in addition to signal fidelity in the subspace basis optimization objective.We approximate the compressed domain CRB with a computationally efficient alternative, which yields a sufficiently optimal solution in practice.The proposed CRB-SVD basis is a drop-in replacement for the traditional SVD basis with improved representational compactness, promoting improved parameter accuracy and precision at smaller subspace sizes and offering computational speedups and memory savings for subspace reconstruction tasks.

Fig. 1 .
Fig. 1.(a) Distinguishing one model parameter from another depends on the components of its signal derivative, j i , that are orthogonal to the span of all other signal derivatives, J i[44].(b) Depiction of a representative signal s for the inversion recovery MRF-FISP sequence[26] (assuming the scaling M 0 = 1), its derivative with respect to the longitudinal relaxation rate ∂s/∂R 1 , and the corresponding orthogonalized derivative (∂s/∂R 1 ) ⊥ , where the components parallel to the derivatives with respect to all other model parameters were removed (i.e., M 0 and R 2 ).Note we plot the absolute value and (∂s/∂R 1 ) ⊥ is scaled up to have a unit length, which emphasizes that the first segment is most important for encoding R 1 .

Fig. 2 .
Fig. 2. (a) Signal energy loss vs λ normalized by the total signal energy in the qMT test dataset.(b) Cram ér-Rao bound (CRB) loss calculated using the approximate compressed CRB (∆Bac; Eq. (6)) and the exact compressed CRB (∆Bec; Eq. (8)), averaged over all fingerprints and orthogonalized derivatives in the test dataset.(c) CRB loss for m s 0 and Rx individually at Nc = 15.(d) Average ratio of approximate compressed CRB to exact compressed CRB over the test dataset (Eq.(13)).The λ values shaded in gray offer significant CRB improvements at a small cost to the signal fidelity for all Nc.

Fig. 3 .
Fig. 3. Representative scatter plots of the approximate compressed CRB of the semi-solid spin pool fraction (Bac(m s 0 )) vs the exact compressed CRB (Bec(m s 0 )) for different numbers of coefficients (Nc) and CRB-loss weightings (λ), where each dot represents a fingerprint in the qMT test set.Note each subplot's axes are limited to the same range for improved visualization.Linear regressions across all data points in each subplot (including those not plotted within the axis limits) are shown in purple in reference to the identity line in red.Correlations between Bac and Bec are improved with increasing Nc and λ.

Fig. 4 .
Fig. 4. Bias (a-b) and standard deviation (c-d) of non-linear least squares (NLLS)-based m s 0 and Rx estimates of a typical white matter fingerprint derived from the traditional SVD (λ = 0; blue) and CRB-SVD (λ = 0.5; red) bases as a function of the number of coefficients Nc.Both metrics are normalized by the ground-truth m s 0 and Rx values, and the reference lines in (c-d) indicate the uncompressed CRB (Eq.(2)) and exact compressed CRB (Eq.(4)).Lower bias and variance more closely resembling minimum variance unbiased estimation is observed with the CRB-SVD basis, particularly for Nc < 18, where the traditional SVD basis performs especially poorly.

Fig. 5 .
Fig. 5. Comparison of non-linear least squares fits of the semi-solid spin pool's fractional size m s 0 (a) and magnetization exchange rate Rx (b) extracted from the traditional SVD (λ = 0) and proposed CRB-SVD basis (λ = 0.5) for different numbers of coefficients (Nc), where Nc = 30 serves as a gold standard.The CRB-SVD basis improves the precision of the m s 0 (magnifications) and the accuracy of the Rx maps for small Nc.

Fig. 6 .
Fig. 6.Comparison of R 1 = 1/T 1 maps extracted from the MRF-FISP data that were reconstructed with a traditional SVD basis (λ = 0) and the proposed CRB-SVD basis (λ = 0.3) for different subspace sizes (Nc), where Nc = 10 serves as a gold standard.The CRB-SVD basis offers improved parameter accuracy at Nc = 3, the number of parameters in the biophysical model, and is non-inferior for Nc ≥ 4.